%  function [sig_active_cells mfr] = analyze_mfr_fmri( start_time_msecs, end_time_msecs,significance_level, save_seconds)
%
% Reads from spikes*.dat files in the current directory of the form: 
% 
% 0 0
% 250 1
% 500 2
% 750 3
% 1000 0
% 1250 1
% 1500 2
% 1750 3
% 2000 0
% 2250 1
% 2500 2
% 2750 3
% ...
% 
% 
% mfr is defined as follows: mfr is the mean firing rate of all cells in
% the network over the whole time period.

function [mfr] = analyze_mfr_fmri( start_time_msecs, end_time_msecs,significance_level, save_seconds)

%global samples cats
global mfr rows cols
msecs_per_sec = 1000;

if ~exist('significance_level')
    significance_level=0.05;
end
if ~exist('save_seconds')
    save_seconds = 1;
end


read_groups; % Assumes groups.dat has the neural area data

num_neurons = max(neuron_id(:,2))+1

num_seconds = (end_time_msecs - start_time_msecs + 1)/1000

spike_counts=zeros(num_neurons, num_seconds);


filename = '';

second = 1;
for t = start_time_msecs:save_seconds*msecs_per_sec:end_time_msecs
    
    % find the file we need to open
    file = ceil( (t+1)/(save_seconds*1000));

    filename = ['spikes' num2str( (file)*save_seconds*1000, 15 ) '.dat']
    sdat = load(filename);

    spike_counts(:,second) = hist( sdat(:,2), 0:(num_neurons-1))'; % Make a bin for each neuron.

    second = second + 1;
end

mfr = sum(spike_counts,2)*(1/num_seconds);
stdev = std(spike_counts,0,2);
dof = num_seconds-1;
%sig_active_cells = zscore(spike_counts')'>tinv(1-significance_level,dof);
%for i = 1:1000,[h(i) p(i)]=ttest2(a(1:1000),a(0+i),.001);end

%[junk significant]=ranksum( mfr( i ).spikes(n, : ) , mfr( j ).spikes(n, : ) ,significance_level);



